import matplotlib
matplotlib.use('PDF')
from scipy.optimize import leastsq
from optimizationFunctions.residual_enthalpy_data import residual_enthalpy_data
from optimizationFunctions.residual_pXT_data import residual_pXT_data
from optimizationFunctions.calc_enthalpy import calc_enthalpy
from optimizationFunctions.residual_enthalpy_and_PVT import residual_enthalpy_and_PVT
#from Calc_B_from_B12 import Calc_B_from_B12
from pylab import *
from numpy import *
import xlrd
import numpy

MPa_to_kPa=1000.0
Bars_to_kPa=100.0
   
p=[ 3.93511381e-01,   5.91798638e-01,   2.17340160e+00,   6.58736091e-01,
   1.49527903e+00,   1.80554325e+00,   1.00000975e+00,
   2.15932933e+00,  1.70370289e+00,   9.99993200e-01,   1.71161665e+00,   
   2.33257272e-05, 1.06367793e+01,   1.35999949e+00,
   5.81075049e-01,   2.00846269e+01, 2.80172838e+00,   2.67567649e+00,   
   1.85430123e-04,   3.11521103e-01, 1.92272528e+00,   
   1.48352080e+00,  -2.27589359e-01,   1.97388639e+00, 1.40630453e+00,
   2.01963999e-02,   9.61708221e-01,   5.09904925e+00, 
   1.52479618e-01, 6.79051907e-01,   9.66836715e+00,   2.37141683e+00]
   

p=asarray(p)

#molecular weight of h2 and h2o
h2_M_amu=2.01594
h2o_M_amu=float(8.314472/0.46151805)

fname='H2-H2O_EOS_data.xls'
#Read Experimental Data From Excel xls
wb=xlrd.open_workbook(fname)
pvtSheet=wb.sheet_by_index(0) #pVT from Seward and Franck 1981
enthalpySheet=wb.sheet_by_index(1) # enthalpy data from Lancaster and Wormald, 1990
bcmeasSheet=wb.sheet_by_index(3) # B, C virial coef from Seward 2000
synFrombcSheet=wb.sheet_by_index(2) # synthetic pVT data set taken from Seward 2000
btheorySheet=wb.sheet_by_index(4) # B12 data from Hodges 2004
 
Data_Enthalpy_x_h2=asarray(enthalpySheet.col_values(0)[1:len(enthalpySheet.col_values(0))])
Data_Enthalpy_x_h2o=1.0-Data_Enthalpy_x_h2
Data_Enthalpy_T=asarray(enthalpySheet.col_values(1)[1:len(enthalpySheet.col_values(1))])
Data_Enthalpy_PMPa=asarray(enthalpySheet.col_values(2)[1:len(enthalpySheet.col_values(2))])
Data_Enthalpy_H=asarray(enthalpySheet.col_values(3)[1:len(enthalpySheet.col_values(3))])
Data_Enthalpy_PkPa=asarray(Data_Enthalpy_PMPa)*MPa_to_kPa

Data_pXT_T=asarray(pvtSheet.col_values(0)[1:len(pvtSheet.col_values(0))])+273.15 #Celcius-->Kelvin
Data_pXT_P=asarray(pvtSheet.col_values(1)[1:len(pvtSheet.col_values(1))])*Bars_to_kPa # Bars->kPa
Data_pXT_x_h2=asarray(pvtSheet.col_values(2)[1:len(pvtSheet.col_values(2))])

Data_pXT_x_h2o=1.0-Data_pXT_x_h2 
Data_pXT_Mmix=((Data_pXT_x_h2o*h2o_M_amu)+(Data_pXT_x_h2*h2_M_amu))
Data_pXT_Vol=asarray(pvtSheet.col_values(3)[1:len(pvtSheet.col_values(3))])*(1.0/Data_pXT_Mmix)*(1.0/1.0e6)*1000.0 # cm**3/mol ->m**3/kg
                                                                                                          #(cm**3/mol)(1 m**3/1e6 cm**3)(1 mol/Mmix(g))(1000 g/ 1 kg)


Data_pXT_syn_xh2=asarray(synFrombcSheet.col_values(0)[1:len(synFrombcSheet.col_values(0))])/100.
Data_pXT_syn_xh2o=1-Data_pXT_syn_xh2

Data_pXT_syn_T=asarray(synFrombcSheet.col_values(3)[1:len(synFrombcSheet.col_values(3))])+273.15
Data_pXT_syn_P=asarray(synFrombcSheet.col_values(4)[1:len(synFrombcSheet.col_values(4))])*Bars_to_kPa # Bars->kPa

Data_pXT_Mmix=((Data_pXT_syn_xh2o*h2o_M_amu)+(Data_pXT_syn_xh2*h2_M_amu))
Data_pXT_syn_Vol=asarray(synFrombcSheet.col_values(2)[1:len(synFrombcSheet.col_values(2))])*(1.0/Data_pXT_Mmix)*(1.0/1.0e6)*1000.0 # cm**3/mol ->m**3/kg
                                                                                                          #(cm**3/mol)(1 m**3/1e6 cm**3)(1 mol/Mmix(g))(1000 g/ 1 kg)
Data_pXT_T=numpy.append(Data_pXT_T,Data_pXT_syn_T)
Data_pXT_P=numpy.append(Data_pXT_P,Data_pXT_syn_P)
Data_pXT_Vol=numpy.append(Data_pXT_Vol,Data_pXT_syn_Vol)
Data_pXT_x_h2=numpy.append(Data_pXT_x_h2,Data_pXT_syn_xh2)
Data_pXT_x_h2o=numpy.append(Data_pXT_x_h2o,Data_pXT_syn_xh2o)

Data5050_Enthalpy_T=Data_Enthalpy_T[0:61]
Data5050_Enthalpy_P=Data_Enthalpy_PkPa[0:61]/Bars_to_kPa
Data5050_Data_Enthalpy_H=Data_Enthalpy_H[0:61]

sorted_index=Data5050_Enthalpy_T.argsort(axis=0)
Enthalpy_T_sorted=Data5050_Enthalpy_T[sorted_index]

Enthalpy_P_sorted=Data5050_Enthalpy_P[sorted_index]
Enthalpy_H_sorted=Data5050_Data_Enthalpy_H[sorted_index]
##
#PVT Laboratory Data taken by Karpowicz and Steffes
##
myDataT=[376.204932,377.726416,446.697607]
myDataP=[19.367289432236020*Bars_to_kPa,75.059974199647826*Bars_to_kPa,87.860145913789353*Bars_to_kPa]
myDataXH2=[18.663980319129138/(18.663980319129138+1.110081400053644),74.964042318807220/(74.964042318807220+1.110081400053644),74.964042318807220/(74.964042318807220+1.110081400053644)]                  

myDataDensity=[] #g/L equivalent to kg/m^3
myDataDensity.append((18.663980319129138*h2_M_amu+1.126616830674388*h2o_M_amu)/32.328408091385590)
myDataDensity.append((74.965164072252165*h2_M_amu+1.126616830674388*h2o_M_amu)/32.337938041626330)
myDataDensity.append((74.965164072252165*h2_M_amu+1.126616830674388*h2o_M_amu)/32.413598694129391)
##

Data_pXT_density=1.0/Data_pXT_Vol #kg/m**3                      
Data_pXT_T=numpy.append(Data_pXT_T,numpy.asarray(myDataT))
Data_pXT_P=numpy.append(Data_pXT_P,numpy.asarray(myDataP))      
Data_pXT_x_h2=numpy.append(Data_pXT_x_h2,numpy.asarray(myDataXH2)) #g/L equivalent to kg/m^3
Data_pXT_x_h2o=1.0-Data_pXT_x_h2
Data_pXT_density=numpy.append(Data_pXT_density,numpy.asarray(myDataDensity))


CC_per_mol_to_Liters_per_mol=0.001



Data_B_T_vals=273.15+asarray(bcmeasSheet.col_values(0)[1:len(bcmeasSheet.col_values(0))]) #in C converted to K
Data_B_x_h2_vals=asarray(bcmeasSheet.col_values(1)[1:len(bcmeasSheet.col_values(1))])
Data_B_x_h2o_vals=1.0-Data_B_x_h2_vals
Data_B_B_vals=CC_per_mol_to_Liters_per_mol*asarray(bcmeasSheet.col_values(2)[1:len(bcmeasSheet.col_values(2))])
Data_C_vals=CC_per_mol_to_Liters_per_mol*CC_per_mol_to_Liters_per_mol*asarray(bcmeasSheet.col_values(3)[1:len(bcmeasSheet.col_values(3))])                       

#Sythetic Bmix data using B12 from Hodges, 2004
Data_B_theory_B=CC_per_mol_to_Liters_per_mol*asarray(btheorySheet.col_values(4)[1:len(btheorySheet.col_values(2))])
Data_B_theory_T=asarray(btheorySheet.col_values(0)[1:len(btheorySheet.col_values(0))])                                                  
Data_B_theory_xh2=asarray(btheorySheet.col_values(5)[1:len(btheorySheet.col_values(5))])
Data_B_theory_xh2o=1.-Data_B_theory_xh2                           
opt=1

[optimized_p,message]=leastsq(residual_enthalpy_and_PVT,p[:],args=(Data_Enthalpy_H,
                                                                   Data_Enthalpy_PkPa,
                                                                   Data_Enthalpy_T,
                                                                   Data_Enthalpy_x_h2,
                                                                   Data_Enthalpy_x_h2o,
                                                                   Data_pXT_P,
                                                                   Data_pXT_T,
                                                                   Data_pXT_density,
                                                                   Data_pXT_x_h2,
                                                                   Data_pXT_x_h2o,
                                                                   Data_B_B_vals,
                                                                   Data_B_T_vals,
                                                                   Data_B_x_h2_vals,
                                                                   Data_B_x_h2o_vals,
                                                                   Data_C_vals,
                                                                   Data_B_theory_B,
                                                                   Data_B_theory_T,
                                                                   Data_B_theory_xh2,
                                                                   Data_B_theory_xh2o,opt),maxfev=10000000)

#optimized_p=p
opt=0
#print 'final optimized',optimized_p
print 'BetaT', abs(p[0])
print 'BetaV',abs(p[1])
print 'GammaT',abs(p[2])         
print 'GammaV',abs(p[3])

print 'Nk poly',abs(p[4]),abs(p[11]),abs(p[18]),abs(p[25])
print 'tk Poly',p[5],p[12],p[19],p[26]
print 'dk poly',p[6],p[13],p[20],p[27]

print 'Nl exp',abs(p[7]),abs(p[14]),abs(p[21]),abs(p[28])
print 'tl exp',p[8],p[15],p[22],p[29]
print 'dl exp',p[9],p[16],p[23],p[30]
print 'ck exp',p[10],p[17],p[24],p[31]



p=optimized_p

x_h2_598_2=asarray(arange(0.01,1.0,0.01))
x_h2_698_2=asarray(arange(0.01,1.0,0.01))


P448=asarray(arange(0.1,1.3,0.1)*MPa_to_kPa)
P473=asarray(arange(0.1,2,0.1)*MPa_to_kPa)
P498=asarray(arange(0.1,3,0.1)*MPa_to_kPa)
P523=asarray(arange(0.1,4,0.1)*MPa_to_kPa)
P548=asarray(arange(0.1,6,0.1)*MPa_to_kPa)
P573=asarray(arange(0.1,8.5,0.1)*MPa_to_kPa)
P598=asarray(arange(0.1,11,0.1)*MPa_to_kPa)
P648=asarray(arange(0.1,13,0.1)*MPa_to_kPa)
P698=asarray(arange(0.1,13,0.1)*MPa_to_kPa)

P_598_2=10.51*MPa_to_kPa*ones(shape(x_h2_598_2))
P_698_2=11.13*MPa_to_kPa*ones(shape(x_h2_698_2))

T448=448.2*ones(shape(P448))
T473=473.2*ones(shape(P473))
T498=498.2*ones(shape(P498))
T523=523.2*ones(shape(P523))
T548=548.2*ones(shape(P548))
T573=573.2*ones(shape(P573))
T598=598.2*ones(shape(P598))
T648=648.2*ones(shape(P648))
T698=698.2*ones(shape(P698))

T_598_2=598.2*ones(shape(P_598_2))
T_698_2=698.2*ones(shape(P_698_2))

x_h2_448=0.5*ones(shape(P448))
x_h2_473=0.5*ones(shape(P473))
x_h2_498=0.5*ones(shape(P498))
x_h2_523=0.5*ones(shape(P523))
x_h2_548=0.5*ones(shape(P548))
x_h2_573=0.5*ones(shape(P573))
x_h2_598=0.5*ones(shape(P598))
x_h2_648=0.5*ones(shape(P648))
x_h2_698=0.5*ones(shape(P698))

x_h2o_698_2=1.0-x_h2_698_2
x_h2o_598_2=1.0-x_h2_598_2

Curve_448=calc_enthalpy(p,P448,T448,x_h2_448,x_h2_448)
Curve_473=calc_enthalpy(p,P473,T473,x_h2_473,x_h2_473)
Curve_498=calc_enthalpy(p,P498,T498,x_h2_498,x_h2_498)
Curve_523=calc_enthalpy(p,P523,T523,x_h2_523,x_h2_523)
Curve_548=calc_enthalpy(p,P548,T548,x_h2_548,x_h2_548)
Curve_573=calc_enthalpy(p,P573,T573,x_h2_573,x_h2_573)
Curve_598=calc_enthalpy(p,P598,T598,x_h2_598,x_h2_598)
Curve_648=calc_enthalpy(p,P648,T648,x_h2_648,x_h2_648)
Curve_698=calc_enthalpy(p,P698,T698,x_h2_698,x_h2_698)

Curve_698_2=calc_enthalpy(p,P_698_2,T_698_2,x_h2_698_2,x_h2o_698_2)
Curve_598_2=calc_enthalpy(p,P_598_2,T_598_2,x_h2_598_2,x_h2o_598_2)

figure(1)
plot(P448/Bars_to_kPa,Curve_448,'r--',label='448 K')
plot(P473/Bars_to_kPa,Curve_473,'g--',label='473 K')
plot(P498/Bars_to_kPa,Curve_498,'b',label='498 K')
plot(P523/Bars_to_kPa,Curve_523,'k',label='523 K')

plot(P548/Bars_to_kPa,Curve_548,'c', label='548 K')
plot(P573/Bars_to_kPa,Curve_573,'m', label='573 K')
plot(P598/Bars_to_kPa,Curve_598,'y', label='598 K')
plot(P648/Bars_to_kPa,Curve_648,'r', label='648 K')
plot(P698/Bars_to_kPa,Curve_698,'g', label='698 K')

errorbar(Enthalpy_P_sorted[0:5],Enthalpy_H_sorted[0:5],yerr=0.02*Enthalpy_H_sorted[0:5],ecolor='r',fmt=None)
errorbar(Enthalpy_P_sorted[5:10],Enthalpy_H_sorted[5:10],yerr=0.02*Enthalpy_H_sorted[5:10],ecolor='g',fmt=None)
errorbar(Enthalpy_P_sorted[10:15],Enthalpy_H_sorted[10:15],yerr=0.02*Enthalpy_H_sorted[10:15],ecolor='b',fmt=None)
errorbar(Enthalpy_P_sorted[15:20],Enthalpy_H_sorted[15:20],yerr=0.02*Enthalpy_H_sorted[15:20],ecolor='k',fmt=None)

errorbar(Enthalpy_P_sorted[20:25],Enthalpy_H_sorted[20:25],yerr=0.02*Enthalpy_H_sorted[20:25],ecolor='c',fmt=None)
errorbar(Enthalpy_P_sorted[25:30],Enthalpy_H_sorted[25:30],yerr=0.02*Enthalpy_H_sorted[25:30],ecolor='m',fmt=None)
errorbar(Enthalpy_P_sorted[30:40],Enthalpy_H_sorted[30:40],yerr=0.02*Enthalpy_H_sorted[30:40],ecolor='y',fmt=None)
errorbar(Enthalpy_P_sorted[40:51],Enthalpy_H_sorted[40:51],yerr=0.02*Enthalpy_H_sorted[40:51],ecolor='r',fmt=None)
errorbar(Enthalpy_P_sorted[51:61],Enthalpy_H_sorted[51:61],yerr=0.02*Enthalpy_H_sorted[51:61],ecolor='g',fmt=None)
legend(loc='upper left')
xlabel('Pressure (bars)')
ylabel('Excess Enthalpy (J/Mol)')
title('Excess Enthalpy 50/50 H$_2$/H$_2$O')
savefig('Excess_Enthalpy_5050.pdf',papertype='letter')
start=len(Data_Enthalpy_x_h2)-10
end=len(Data_Enthalpy_x_h2)-5

figure(2)
plot(1-x_h2o_598_2,Curve_598_2,'b',label='598.2 K, 10.51 MPa')
plot(1-x_h2o_698_2,Curve_698_2,'r',label='698.2 K, 11.13 MPa' )
errorbar(Data_Enthalpy_x_h2[start:end],Data_Enthalpy_H[start:end],yerr=0.02*Data_Enthalpy_H[start:end],ecolor='b',fmt=None)
start=len(Data_Enthalpy_x_h2)-5
end=len(Data_Enthalpy_x_h2)
errorbar(Data_Enthalpy_x_h2[start:end],Data_Enthalpy_H[start:end],yerr=0.02*Data_Enthalpy_H[start:end],ecolor='r',fmt=None)

legend(loc='upper right')
xlabel('Mole Fraction X$_{H_2}$')
ylabel('Excess Enthalpy (J/mol)')
title('Excess Enthalpy for Increasing X$_{H_2}$')
savefig('Excess_Enthalpy_mole_frac.pdf',papertype='letter')
residuals=residual_enthalpy_and_PVT(p[:],Data_Enthalpy_H,Data_Enthalpy_PkPa,Data_Enthalpy_T,Data_Enthalpy_x_h2,Data_Enthalpy_x_h2o,Data_pXT_P,Data_pXT_T,Data_pXT_density,Data_pXT_x_h2,Data_pXT_x_h2o,Data_B_B_vals,Data_B_T_vals,Data_B_x_h2_vals,Data_B_x_h2o_vals,Data_C_vals,
                                    Data_B_theory_B,
                                    Data_B_theory_T,
                                    Data_B_theory_xh2,
                                    Data_B_theory_xh2o,
                                    opt)

figure(3)
#residuals[0:3]=residuals[0:3]
plot(Data_pXT_T,residuals[0:len(Data_pXT_P)],'kx')
xlabel('Temperature (K)')
ylabel('Error in Pressure (%)')
#ylim(-100,100)
title('Residual Pressure from EOS Fit')
savefig('Residual_Pressure.pdf',papertype='letter')

begin=len(Data_pXT_P)+len(Data_Enthalpy_H)
end=begin+len(Data_B_T_vals)

figure(4)
plot(Data_B_T_vals,residuals[begin:end],'kx')
xlabel('Temperature (K)')
ylabel('Second Virial Coefficient Error (%)')
title('Residual Second Virial Coefficient from EOS Fit')
savefig('Residual_B.pdf',papertype='letter')
print 'final optimized',optimized_p
figure(9)
plot(Data_B_x_h2_vals,residuals[begin:end],'kx')
xlabel('Mole Fraction Hydrogen')
ylabel('Second Virial Coefficient Error (%)')
title('Residual Second Virial Coefficient from EOS Fit')
savefig('Residual_B_mole.pdf',papertype='letter')

begin=end
end=begin+len(Data_B_T_vals)
figure(5)
plot(Data_B_T_vals,residuals[begin:end],'kx')
xlabel('Temperature (K)')
ylabel('Third Virial Coefficient Error (%)')
title('Residual Third Virial Coefficient from EOS Fit')
savefig('Residual_C.pdf',papertype='letter')

begin=end
figure(6)
end=begin+len(Data_B_theory_B)
plot(Data_B_theory_T,residuals[begin:end],'kx')
xlabel('Temperature (K)')
ylabel('Second Virial Coefficient Error (%)')
#ylim([-200,200])
title('Error Second Virial Coefficient (Synthetic Data from Hodges, 2004)')
savefig('Residual_B_theory.pdf',papertype='letter')
figure(7)
plot(Data_B_theory_xh2,residuals[begin:end],'kx')
xlabel('Mole Fraction H$_2$')
ylabel('Second Virial Coefficient Error (%)')
#ylim([-200,200])
title('Error Second Virial Coefficient (Synthetic Data from Hodges, 2004)')
savefig('Residual_B_theory_mole_h2.pdf',papertype='letter')
print 'BetaT', abs(p[0])
print 'BetaV',abs(p[1])
print 'GammaT',abs(p[2])         
print 'GammaV',abs(p[3])

print 'Nk poly',abs(p[4]),abs(p[11]),abs(p[18]),abs(p[25])
print 'tk Poly',p[5],p[12],p[19],p[26]
print 'dk poly',p[6],p[13],p[20],p[27]

print 'Nl exp',abs(p[7]),abs(p[14]),abs(p[21]),abs(p[28])
print 'tl exp',p[8],p[15],p[22],p[29]
print 'dl exp',p[9],p[16],p[23],p[30]
print 'ck exp',p[10],p[17],p[24],p[31]
